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Abstract 



We consider the stationary Hamilton- Jacobi equation 



b i3 {x)u Xz u X] = [f(x)} 2 



where b can vanish at some points, the right-hand side / is strictly 
positive and is allowed to be discontinuous. More precisely, we consider 
special class of discontinuities for which the notion of viscosity solution 
is well-suited. We propose a semi-Lagrangian scheme for the numerical 
approximation of the viscosity solution in the sense of Ishii and we 
study its properties. We also prove an a-priori error estimate for the 
scheme in L 1 (fi). The last section contains some applications to control 
and image processing problems. 

Hamilton- Jacobi equation, discontinuous Hamiltonian, viscosity 
solutions, semi- Lagrangian schemes, a-priori error estimates. 
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1 Introduction 

In this paper we study the following boundary value problem. Let Q C l w 
be an open bounded domain with a Lipschitz boundary dQ, we consider the 
Dirichlet problem 



where / and g are given functions. We focus to the fact that the function / 
is Borel measurable, possibly discontinuous. 
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In the most classical case, the matrix (fry) is the identity matrix and / 
is a positive function, so the partial differential equation in ([I]) reduces to 

\Du(x)\ = f(x). (2) 

which is the classical form of an eikonal equation. 

This equation arises in the study of many problems, e.g. in geometrical 
optics, computer vision, control theory and robotic navigation. In geometri- 
cal optics, to describe the propagation of light the eikonal equation appears 
in the form 

N 

bij(x)u Xi u Xj (x) = f{x) (3) 

»J=1 

where b = oo l and / has the meaning of the refraction index of the media 
where the light rays are passing. Typically, the refraction law applies across 
surfaces of discontinuity of /. 

Another example is offered by a classical problem in computer vision, the 
Shape-from-Shading model. In this classical inverse problem we come up 
with the equation 

Vl + \Du(x)\* = j^ (4) 

if we assume that the light source is vertical and at infinity (so all the rays are 
parallel) and the object to reconstruct is the graph of the unknown function 
u. In this particular case the brightness I(x) G (0,1], i.e. the intensity of 
light reflected by the object, can be discontinuous when the object has edges 
as we will see in more details at the end of this paper. Using classical tools 
of convex analysis, both equations above can be rewritten in the form Q. 

Another motivation to deal with discontinuous hamiltonians comes di- 
rectly from control theory. In this framework discontinuous functions can be 
used to represent targets (for example using / as a characteristic function) 
and/or state constraints (using / as an indicator function) [7J. Clearly, the 
well-posedness of ([I]) in the case of continuous / follows from the theory 
of viscosity solutions for HJ equations, the interested reader can find the 
details in [3] and [2] where there are summarized the well-known results 
introduced by Crandall, Lions, Ishii and other authors. It is interesting to 
point out that, when the hamiltonian is discontinuous, the knowledge of / 
at every point will not guarantee the well-posedness of the problem even 
in the framework of viscosity solutions. In fact, for equation ([!]) it can be 
easily observed that, even when / is defined point wise and has appropriate 
discontinuities, the value function for the corresponding control problem will 
not satisfy the equation in the viscosity sense. In order to define viscosity 
solutions for this case, we use appropriate semicontinuous envelopes of /, 
following some ideas introduced by Ishii in |18| . 

The notion of viscosity solution in the case of discontinuous Hamiltonian 
was proposed by Ishii in |18] where he presents some existence and regularity 
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results. Other results of well-posedness of Hamilton-Jacobi equations in 
presence of discontinuous coefficients are been presented by various authors 
in several works (see HS1 HI E] ) and in the specific case, present in many 
applications, of the eikonal equation (32J, [23] . 

Our primary goal is to prove convergence for a semi-Lagrangian scheme 
which has shown to be rather effective in the approximation of Hamilton- 
Jacobi equations. The results which have been proved for this type of 
schemes work for convex and non convex hamiltonians but use the uniform 
continuity of the hamitonian. Moreover, the typical convergence result is 
given for the L°° norm which is rather natural when dealing with classical 
viscosity solutions (see e.g. the result by Crandall and Lions [9], Barles and 
Souganidis [S] and the monograph by Falcone and Ferretti [E]). For clas- 
sical viscosity solutions, at our knowledge, the only two convergence results 
in L 1 (0) has been proved by Lin and Tadmor [2"T] for a central finite 
difference scheme and by Bokanowwsky et al. [7j in dimension one. We have 
also to mention the level set approach for discontinuous solutions proposed 
by Tsai et al. |31j . Although classical schemes tailored for the the approxi- 
mation of regular cases with convex hamiltonians can give reasonable results 
also for some discontinuous hamiltonians, it would be interesting to have a 
theoretical results in this situation. Deckelnick and Elliott [13] have stud- 
ied a problem where the solution is still Lipschitz continuous although the 
hamiltonian is discontinuous. In particular, they have proposed a finite dif- 
ference scheme for the approximation of ([2]) and their scheme is very similar 
to a finite difference schemes usually applied for regular hamiltonians. Their 
contribution is interesting because they prove an a-priori error estimate in 

Although our work has been also inspired by their results, we use differ- 
ent techniques and our analysis is devoted to a scheme of semi-Lagrangian 
type. The benefits of a semi-Lagrangian scheme with respect to a finite 
differences scheme are a better ability to follow the informations driven by 
the characteristics, the fact that one can use a larger time-step in evolutive 
problems still having stability and the fact that SL-schemes do not require 
a structured grid. This peculiarities give us a faster and more accurate ap- 
proximation in many cases as it has been reported in the literature (see e.g. 
|14| [TT] or appendix A of |2]). It is also important to note that we prove 
an a-priori error estimate which improves the result in |13j because we con- 
sider a more general case ([I]) where also discontinuous viscosity solutions 
are accepted. 

This paper is organized as follows. 
In Section [2] we recall some definitions and theoretical results available for 
discontinuous hamiltonian. Section [3] is devoted to the presentation of the 
scheme and to the proof of some properties which will be used in the proof 
of convergence. In Section [4] we prove convergence and establish an a-priori 
error estimate giving the rate of convergence in the L l (£l) norm. Finally, 
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in Section [5] we present our numerical experiments dealing with control and 
image processing problems. 

2 The model problem and previous results 

We present, for readers convenience, some results of well-posedness mainly 
taken from a work of Soravia [28]. We also introduce our assumptions, which 
are summarized below. 
The boundary data 

g : dCl — > [0, +oo[ is continuous, (5) 

the matrix of the coefficients satisfies 

(kj) = (a ik ) ■ (a* kj ) (6) 

where i,j = 1,...,N and k = 1,...,M and (M < N). Then (bij) is a 
symmetric, positive semidefinite and possibly degenerate matrix, 

cr(-) = {(Jik)i=i,...N; k=i,...M '■ ^ — > is L-Lipschitz continuous. (7) 

moreover, the function / : M. N — > [p, +oo[, p > is Borel measurable and 
possibly discontinuous. 

We can give an optimal control interpretation of ([T]), rewriting the dif- 
ferential operator in the following form 

M 

bij(x)pipj = ^(p ■ a k {x)) 2 = \p ■ a(x)\ 2 , (8) 

k=l 

where the columns of the matrix {<Jik)i,k are the vector fields which will be 
denoted by a k : R N , k = 1, ...M. We define M a = maxj SfcCTj fc. In 
this way the eikonal equation Q becomes, for a = (a 1 , ...a M ) E M. M , the 
following Bellman equation 

max(-Du(x) • V a k a k (x)\ = f(x) (9) 

' a '- 1 1 ti J 

associated to the symmetric controlled dynamics 

M 

y = Y,a k a k (y), y(0) = x, (10) 
k=l 

where the measurable functions a : [0, +oo[— > {a G M. M : \a\ < 1} are the 
controls. We will denote in the sequel by y x {-) '■= Vx(-,cl) the solutions of 



(10). In this system, optimal trajectories are the geodesies associated to the 
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metric defined by the matrix (6f/), and they are straight lines when (bij) 
is the identity matrix. Solution of the equation ^ minimize the following 
functional 

J(x, o(-)) = r f(v(t))dt + g(y(r x )) (11) 



where T x (a(-)) = inf{t : y x (t,a) ^ fi}. 

Let us introduce the concept of discontinuous viscosity solution for ([I]) 
introduced by Ishii in [15] . 

Let / be bounded in £1 and let 

f*(x) = lim inf{/(y) : |y - x\ < r} (12) 

r->0+ 

f*(x) = lim sup{/(y) : \y - x\ < r} (13) 

/* and /* are respectively the lower semicontinuous and the upper semicon- 
tinuous envelope of /. 

Definition 1. A lower (resp. upper) semicontinuous function u : 0, — > 
MU{+oo} (resp. u : £1 ^ M) is a viscosity super- (resp. sub-) solution of the 
equation |7]) if for every if £ C 1 ^), u{x) < +oo, andx G argmin xen (u—(p), 
(resp. x £ argmax xe ^(u — ip) ), we have 

bij(x)(p Xt (x)ip Xj (x) > [f*(x)} 2 , (resp. bij(x)(p Xi (x)ip Xj (x) < [f*{x)f ). 

A function u is a discontinuous viscosity solution of |7p ifu* is a subsolution 
and u* is a supersolution. 

We remind also that the Dirichlet condition is satisfied in the following 
weaker sense 

Definition 2. An upper semicontinuous function u : U — > K, subsolution of 
|I]), satisfies the Dirichlet type boundary condition in the viscosity sense if 
for all ip G C 1 and x G dQ, x G argmax x& ^{u — cp) such that u(x) > g(x), 
then we have 

bij(x)(p Xi <p x . < [f*{x)} 2 . 

Lower semicontinuous functions that satisfy a Dirichlet type boundary con- 
dition are defined accordingly. 

In order to see how easily uniqueness can fail without proper assumptions 
on /, now that we accepted that envelopes of function should be used let us 
consider the ID equation 

\u'(x)\ = f(x), a: €[-2,2], u(-2) = u(2) = 0. (14) 

with the choice f(x) = 2%q, where xq is the characteristic function of the 
rationals. Then one easily checks that both u\ = and U2 = 2 — 2\x\ are 
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viscosity solutions. It is clear, that in general we do not have uniqueness 
of the discontinuous viscosity solution. We add a key assumption on the 
coefficient /. 

Assumption Al. Let us assume that there exist rj > and K > such 
that for every igSJ there is a direction n = n x E S"™ -1 with 

f(y + rd)-f(y)<Kr (15) 

for every y 6 Q, d G S 71 ^ 1 , r > with |y — sg| < n, \d — n\ < r\ and y + rd G f&. 

Under Assumption Al the following comparison theorem holds. This 
result, under some more general hypotheses, is presented in [28J. 

Theorem 1. Let Q be an open domain with Lipschitz boundary. Assume 
§) 3 §j, and @. Let u, v : — > R 6e respectively an upper and a lower- 
semicontinuous function, bounded from below, respectively a subsolution and 
a supersolution of 

bi,j(x)u Xt u Xj = [f(x)] 2 , 

Let us assume that v restricted to d£l is continuous and that u satisfies the 
Dirichlet type boundary condition. Suppose moreover that u or v is Lipschitz 
continuous. Then u < v in £1. 

Prom this result, it follows directly that we have uniqueness of a contin- 
uous solution. 

Corollary 1. Assume @), §j ; and (Al). Let u : O — > R be a continu- 
ous, bounded viscosity solution of the problem Then u is unique in the 
class of discontinuous solutions of the corresponding Dirichlet type problem. 

Example 1 (Soravia |27J). This example shows that discontinuous solutions 
may exists without any contradiction with the previous result. This is due 
to the fact that Corollary [I] does not cover all possible situations. Let us 
consider the Dirichlet problem 

x 2 (u x (x, y)) 2 + (u y (x, y)) 2 = [f(x, y)} 2 ] - 1, l[x] - 1, 1[ 

u(±l,y) = u(a;,±l) = x,ye[-l,l] { ' 

where f(x,y) = 2, for x > 0, and f(x,y) = 1 for x < 0. In this case we 
have that 



b ^={l l)' ^=(0 I) 



therefore the Bellman's equation in this case is 



m&x{-Du(x,y) ■ ai(x,0) T - Du(x,y) ■ a 2 (0, 1) T } = f(x,y). (17) 

|o|<l 
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It is easy to verify that the piecewise continuous function, 

[1 - \y\) x > 0, 1 2/| > 1 + lnx 
u(x,y) = { -2ln(x) x> 0, jy|<l + lnx (18) 

x < 0. 



u{-x,y) 
2 



is a viscosity solution of the problem. We know, as indirect implication of 
Corollary^ that there is no continuous solution. We note that all the class of 
functions with values in x = between 1 — |y| and 2(1 — \y\) are discontinuous 
viscosity solutions. However, we have that all discontinuous solutions have 
u as upper semicontinuous envelope. 

As shown in Example [TJ in general we do not have existence of a con- 
tinuous solution and, in general, we do no have a unique solution. But 
restricting ourselves to a special class of solutions, essentially the case pre- 
sented in the previous example, we can preserve the accuracy of numerical 
approximations and we can also get an error estimate, as we will see in the 
sequel. 

The presence of discontinuities is due to the degeneracy of the coefficient 
a. To handle this case we need some additional hypotheses. In this case, 
however, the assumption will be given on the interface of degeneracy of a. 

From here we will restrict ourselves to the case N = 2. 

Let us denote by £(C) the length of a curve C and assume the existence 
of a regular curve So which splits the domain fi in two non degenerating 
parts. Calling r](x) = (rj 1 (x) , n 2 (x)) the unit normal to So on the point 
x 6 So, we state: 

Assumption A2. There exists a curve So C such that, for the points 
x £ So we have 

r/ 1 (x)ai(x) + ri 2 (x)a 2 (x) = 0; 

moreover 

1. p l (x)o~i(x) + p 2 (x)a2(x) ^ for every (p l ,p 2 ) £ B(0, 1) and x ^ So; 

2. £(S ) < +oo. 

3. Let = f^i Ul^UEo, where, in each subset Vtj there is not degeneracy 
of a, we have flj n dQ / for j = {1, 2}. 

We conclude this section with the following result, which can be derived 
by adapting the classical proof by Ishii [19] : 



Theorem 2. Let Q be an open domain with Lipschitz boundary. Assume 
|5j) 3 |^), |?[), [/5| and assumptions A2. Let u : — > R be a bounded viscosity 
solution of the problem |IJ). It is Lipschitz continuous in every set and 

n 2 . 
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Proof. Take a parameter 5 > 0, and define the set 

S 5 := {x G n\B(x,S) nS / 0}; (19) 



we want to study the regularity of the viscosity solution in the set fii \ £,$ = 

<>;. 

In order to describe our boundary assumptions on fl 1 n S# let us define 
L : Ul x fij -»■ K by 

L(x,y):=inf{ ^ N(f*( , y(t)),j'(t))dt\j G W 1,oo ((0, l),nj) (20) 
•/ 

with 7(0) = x, 7(1) = y} 

where 

iV(r,C) :=sup|-(C,p)|max|-p-^aV fc (x) =r||. (21) 



fc=i 



Then we extend the boundary condition to n S5 in the following way: 
g(x) = inf { 5 (y) + L(x,y)} x £ n Sj (22) 

We can claim now, that there exists a viscosity solution G C 0,1 (n*) of 
([!]) with the Dirichlet conditions introduced above. This is proved in Ishii 

We do the same on the set Q2, getting the function u s 2 G C 0,1 (r2 2 )- Now 
the class of functions 

7=r<5 

x G S2 1 



(23) 

I X G ii 2 



in a viscosity solution of Q in fi* U For the arbitrariness of 5 and 
defining u on the discontinuity as said previously we get the thesis. □ 

Which value the solution can assume in So? As shown in Example 2 
and in accordance with the definition of discontinuous viscosity solutions, 
we can choose for x G So every value between -u* and u*. 

We can observe that in this class we can also include a easier case. If 
instead the function a we consider a coefficient c(x) : O — > 1R where c(x) > 
for all x G f2 but which can vanish in some points. In particular, in this case 
we will define So := {x G 0|c(x) = 0} and the previous hypothesis on the 
nature of So reduces to 

£(S ) < +00 and % n dtl / for j = {1, 2}. 
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3 The semi-Lagrangian approximation scheme and 
its properties 

We construct a semi-Lagrangian approximation scheme for the equation ([I]) 
following the approach p3] . 

Introducing the Kruzkov's change of variable, v(x) = 1 - e~ u ^ and 
using Q and Q the problem ([!]) becomes 

\Dv(x)-a(x)\ = f(x)(l-v(x)) xeQ . , 

v(x) = 1 - e-f^ x£dn { ' 

to come back to the original unknown u we can use the inverse transform, 
i.e. u(x) = ln(l — v{x)). 

Let us to observe that since u{x) > 0, we have < v(x) < 1. We can 
write the previous equation in the equivalent way 

{ v{x) + j^\Dv{x) ■ a{x)\ = I xen 

I v (x) = 1 - e~ 9 ^ xedQ [ ' 



We want to build a discrete approximation of (25 ). We pass to the Bellman's 
equation in the following form 

sup J Dfc^o-fc(^) .Dv(x)\ = l-v(x). (26) 



a(EB(0,l) I f( x 

We observe that, in this formulation, it exists a clear interpretation of this 
equation as the value function of an optimization problem of constant run- 
ning cost and discount factor equal to one and the modulus of the velocity 
in the direction a of the dynamics equal to a jffi ■ 



We discretize the left-hand side term of (26) as a directional derivative 



and we arrive to the following discrete problem: 

v h(%) = jh ae f^ 1} { y h " jfc) E fc a k a k (x)^ } + ^ x G fi 

Vh (x) = 1 - e~3( x ) xedfl 

(27) 

where h is a positive real number and we will assume (to simplify the pre- 
sentation) that x — j^y Y^k ak(T k(x) G ^ for every a G B(0, 1). 

We have to remark that for x G Vt and a direction d G dB(0, 1), we always 
can find an a G B(0, 1) such that A- = d and x — a k <Jk(x) G 0, (see 

Figure [T]) because is an open set and we can chose the variable a null to 
remain on x. 



Let introduce a space discretization of (27) yielding a fully discrete 



scheme. We construct a regular triangulation of f2 made by a family of 
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simplices Sj, such that ft = UjSj, denoting x m , m = 1,...,L, the nodes of 
the triangulation, by 

Ax := max diam(Sj) (28) 
j 

the size of the mesh (diam(B) denotes the diameter of the set B) and by G 
the set of the knots of the grid. 
We look for a solution of 

f W(x m ) = a mm^I[W}(x m - j^Za k a k {x m )) + x m G G 

\ W(x m ) = 1 - e~ g{Xm) x m £Gndn 

(29) 

where J[W](x) is a linear interpolation of W on the point x, in the space of 
piecewise linear functions on f2 

W Ax ■- { w : H -> R \w G C(J2) and L>w(:e) = c,- for any x G Sj} . 

Theorem 3. Lei x m - ^ fc a k a k (x m ) G 0, /or ewer?/ x m G G, for any 

a G B(0, l)j so i/iere exists a unique solution W of (29) in W 



Proof. By our assumption, starting from any x m G G we will reach points 
which still belong to Q. So, for every w G W A:r we have 



h X L 



w x r 



f{Xr, 



^2a k a k (x m ) = y^X mj (a)w(xj) 
k J j=i 
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where X m j(a) are the coefficients of the convex combination representing the 
point x m — jfi s a k ak(x m ), and L the number of nodes of G, i.e. 



f{Xm) 



^a fc (T fc (a; m ) = ^ X mj (a)xj 



(30) 



now we observe 



< A m j(a) < 1 and A mj (a) = 1 for any a £ B(0, 1) 



(31) 



Then ( 29 ) is equivalent to the following fixed point problem in finite dimen- 
sion 

W = T{W) 

where the map T : R L R L is defined componentwise as 



(T(W)) m 



1 » / \ T T 7" h 

mm A(a)W + 



mel,...,L (32) 



1 + h aeB(o,i) l + h 

W m = W(x m ) and A(a) is the Lx L matrix of the coefficients X m j satisfying 



(30), (31) for m,j G 1,...,L. 



T is a contraction mapping. In fact, let a be a control giving the mini- 
mum in T(V) m , we have 

[T(W) - T(V)] m < ^ [A(a)(W - V)] m 

< -^—■max\K lj {a)\\\W -V]^ < T ^\\W-V\\ O0 (33) 
1 + a m,j i + n 

Switching the role of W and V we can conclude that 



J_ ||v^ - v\\ 

x - 1 + ^ II Hoc 



(34) 
□ 



3.1 Properties of the scheme 



The solution of (29) has the following crucial proprieties: 
Consistency 



From (29|) we obtain 

W(x m 



\ min \-W(x m )+I[W](x r 
h aeB(o,i) I 



f{Xr, 



S r j a k Q-k{xm)) 



= 1 

(35) 
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We can see the term on the minimum as a first order approximation of the 
directional derivative 

- min {DW -S" a k a k (x) \ + o(h) = 1- W(x m ) (36) 
aeB(0,l)\ ^ j 

using max(-) = — min(— •) we find the consistency, that is of order o(h+Ax). 

Convergence and monotonicity. Since T is a contraction mapping in R , 
the sequence 

W n = T ( W n-l^ (37) 

will converge to W, for any Z E M. N . Moreover, the following estimate holds 
true: n 

||W"-W||oo< (i^) IIWo-^Hoo- (38) 

4 An a-priori estimate in L X (Q) 

In this section we present our main result. As stated in previous section the 
relevance of this result is due to its applicability in the case of discontinuous 
value functions. Using a L x (f2) norm we can extend the convergence result 
also to this class of solutions. 



lions 



Theorem 4. Let assume the hypotheses |5|), |fi]) ; |?p, (i5[ ) and assumpti 
on the set So . Moreover, let < -A- . 
We have that 

\\v(x) - W(x)\\ L i {n) < C^h + C'Ax forallh>0 (39) 
for some positive constant C, C independent from h and Ax. 
Proof. We start introducing the set £ax defined as it follows 
£ Aa , := {x G n\B (x, Ax) n S ^ 0} . 
We observe that 



\v(x) -W(x)\\ L i {n) < / \v{x) -W{x)\dx+ / \v(x) -W(x)\dx 



x)\dx (40) 



<X) / \v(x)-W(x)\dx+ / \v(x)-W( 

j J^Ax 

where f2 := C)jQ,j is the partition of 0, generated from So as stated in the 
definition of the set So- 

Prom the Kruzkov's transform we know that \v(x) — W(x)\ < 1 for all 
x € fl and adding the assumptions on the set Sq we get, for a fixed C > 0, 
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f \v{x) - W(x)\dx < f dx< e(Z )Ax < C'Ax. (41) 

To prove the statement, we need an estimate for the term \v(x)—W(x)\dx 
for every choice of j. With this aim, we remind that, for Theorem [2j both 
v(x) and W{x) are Lipschitz continuous, so we can use a modification of the 
classical argument based on duplication of variables. Something similar can 
be found on [291 [T3] and [28] . 

We are focusing on the problem (29) restricted on the region 0,- := 
Qj \ with some compatible Dirichlet conditions on Qj n d£l. We do not 
have any Dirichlet conditions on dVlj n <9Xaie> so we extend the boundary 
conditions as in (22). Inside the region Qj the solution v(x) is Lipschitz 
continuous by Theorem [2j 

Let us choose a point x G G U Q.j =: Gj such that 

\v(x) - W(x)\ = max \v(x) - W(x)\ (42) 

and assume that v(x) > W(x). The opposite case can be treated similarly. If 
dist(x, dQj) < y/h, implies, from the Dirichlet conditions and the Lipschitz 
continuity of v and W that 



ip(x,y) := v(x)—W(y)—Li — L2^fh\y — x\ 2 , for (x,y) G QjxGj. 



max \v(x) - W{x)\ = v{x) - W{x) < CVh. (43) 

Now, suppose that dist(x,dftj) > \fh and define the auxiliary function 

\x — y — Vht]] 2 

(44) 

Where rj is the inward normal to VI j like stated in previous assumptions. 

It is not hard to check that the boundedness of v, W and the continuity 
of i/j, imply the existence of some (x, y) (depending on h) such that 

^(x,y) > i>{x,y) for all (x,y) G Qj x Gj. (45) 

Since dist(x,dO,j) > \/h, we have that x + y/hrj G flj and therefore 

ip(x, y) > + \/hr), x), (46) 

or equivalently 

v(x)-W(y)-^\x-y-Vhn\ 2 -L 2 Vh\y-x\ 2 > v(x-Vhrj)-W(x). (47) 
v/i 
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(47) implies 



x — y — Vhn\ 2 + Li2*Jh\y — x\ 2 < v(x) — v(x — Vhn) + W{x) — W{y) 



< v(x) — v(y) + [(v(y — W(y)) — {v(x — W{x))} + v(x) — v(x — Vhr/) 
< L v \x — y\ + s/hL v < L v \x — y — Vhr]\ + 2\ThL v 



<^\x-y-Vhr]\ 2 + ^L 2 v + 2VhL v (48) 
2yh 2L\ 



where L v is the Lipschitz constant of v, and therefore we can conclude 
-\x-y-Vh V \ 2 <— 2 Ll + -Ll< 



\y — x\ < 



Li 

2L1L2 v L2 



2 + e 



(49) 
(50) 



for a e > 0, provided L%, L2 are sufficiently large. 

Let us consider now the case (x, x) £ Clj x Gj , so there are not on the 
boundary. 



By (27) we have, for a x £ G j 

^a k a k (x) 



W [x-h- 



W(x) + hW(x) - h 



(51) 



for some a = a(x). This equation is verified a.e. and the point x — 

/(-r)^ ^ from the definition of the admissible choice of a and the 
hypothesis on the discretization steps. Since the map 



x 1 y v(x) 



Wm + L 1 \ X -\-J~ h ^ + L2Vh\y-x\ 2 



has a maximum at x, by (25) we obtain 



and then 



\(x-y-Vhr])-a{x)\ 
Li < f*(x) - U{x)v{x) 



(52) 



(53) 



(54) 



the inequality ip{x, y) > tp [x, y - a k a k {y) ) gives 

\x — y — y/hrjj\ 2 



-W{y)-L x 



2Vh 



f® 

-L 2 Vh\y-x\ 2 > -Wly 



x — hy — vhrj 



2\fh 



LWh 



m 



^a k a k {y) 



y — x — h 



E" gfc(y) 

m 



(55) 
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and then 



>W(y) 



2-Jh 



\x — y — Vhr] 
+ L 2 Vh 



x — y — Vhr] 



E« fcf7 fe(y) 



y-x\ 



y-x 



(56) 



Substituting the left hand side term with (51) and using the fact that for 



every a,b,c G M n we can prove that \a—b\ 2 — \a—b—hc\ 2 = 2h(a—b)-c—h 2 \c\ 2 , 
we get 



2h(x-y-Vh V ). ^ ak{V) -h 2 



+ 



2Vh 



m 

'2 H y-x).^^l-h 2 

f{y) 



(57) 



Now, adding to (54) and using the estimations (49) and fl50 



v(x) - W(y) < Vh + ^ Vh? 

^a k a k (y) *£a k a k (x) 



J2a k a k {y) 



m 



f*{x) 



< \^Vh+ L ^ 



m 

L2\fh(y - x) 

m 



^={x -y- Vhr]) 
^a k a k (x) 



2 + e 



L 2 Vhe 



Finally, choosing e = Vh by the boundedness of / and cr, we obtain 

v(x) - W(y) < CVh 



(58) 



(59) 



where C is a suitable positive constants. Then the inequality vp(x, y) > 
tp(x, x) yields 

v(x) - W{x) < v(x) - W(y) < CVh (60) 



for all x G Q,. 



Finally we consider the case when y G dGj or x G dflj . If y G dGj . the 
Dirichlet conditions imply that v(y) = W(y) and we have 

v(x) — W(x) < v(x — Vhr]) — v{x) + v(y) — v(x) 

< L v (Vh + \x - y\) < L v (2Vh + \x-y- Vhr]\) < CVh. (61) 
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In a similar way we can treat the case x £ dQj. 

To prove the inequality W(x) — v(x) < C\fh it is enough to interchange 
the roles of v and W on the auxiliary function ijj. 



We add this estimation in ( 40 ) , getting the thesis 



\v(x) -W(x)\\ L i < CVh + C'Ax. 



(62) 
□ 



5 Numerical experiments and applications 

In this section we present some results for ([I]) on some test problems coming 
from front propagation, control theory and image processing. In all these 
examples the discontinuity of the coefficients appears in a natural way and 
has an easy interpretation with respect to the model. 

5.1 Test 1: a front propagation problem 

Front propagation problems arise in a lot of different fields of mathematics. 
A typical approach is to use the Hamilton-Jacobi framework to solve them, 
as in [24J. Our first test can be interpreted as a front propagation in a 
discountinuous media. In this model, the level sets of the value function 
have the meaning of the regions with the same time of arrival of the front. 
Let n := (-1, 1) x (0, 2) and / : Q -> M be defined by 



f(xi,x 2 ) 



1 

3/4 
1/2 



xi < 0, 
xi = 
xi > 



(63) 



It is not difficult to see that / satisfies conditions (15). We can verify that 
the function 



u(x 1 ,x 2 ) :-- 



\x 2 , 



xi > 0, 



\x 2 < xi < 0, 



(64) 



X2, 



X\ < 



V3 



X 2 . 



is a viscosity solution of \Du\ = f(x) in the sense of our definition. Moreover, 
we take g := U\qq. We show in the Table [T] and in Figure [2] our results. 

We also show, in Table [2] a comparison between the FD methods pro- 
posed in [13j . They proposed two techniques: in the first there is a regular- 
ization of the Hamiltonian with a viscosity term (DF — reg), in the second 
one (DF — FS), better results are obtained, but numerically there are more 
difficulties; the authors solve them using FastSweeping (see [33J ) as accel- 
eration technique and they archive very good results. Our technique has, in 
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Figure 2: Test 1. 



this test, a performance similar to DF—reg, in our scheme, the interpolation 
operator (in this case bilinear) adds a regularization. We aspect better per- 
formances of our method rather DF techniques on more complicated cases, 
where characteristics are not straight lines. 

5.2 Test 2: a control problem with a discontinuous value 
function 

In this test we present a case where a continuous solution does not exist. 
In this case it is evident that a convergence in uniform norm will not be 
possible. 

We consider the problem shown in the example [TJ As already said, let 
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Ax = h 


1 1 oo 






Ord(Li) 


0.1 


1.734e-l 




8.112e-2 




0.05 


8.039e-2 


1.1095 


3.261e-2 


1.3148 


0.025 


4.359e-2 


0.8830 


1.616e-2 


1.0178 


0.0125 


2.255e-2 


0.9509 


7.985e-3 


1.0271 



Table 1: Test 1: experimental error. 



Ax = h 


our method 


Ord 


DF-reg 


Ord 


DF-FS 


Ord 


0.1 


1.734e-l 




1.243e-l 




5.590e-2 




0.05 


8.039e-2 


1.1095 


7.229e-2 


0.78 


2.795e-2 


1.00 


0.025 


4.359e-2 


0.8830 


4.085e-2 


0.82 


1.397e-2 


1.00 


0.0125 


2.255e-2 


0.9509 


2.266e-2 


0.85 


3.493e-3 


1.00 



Table 2: Test 1: comparison between different numerical methods (uniform 
norm). 



ft := [—1,1] we want to solve 



x 2 (u x (x,y)f + (u y (x,y)r = [f(x,y)} 2 } - 1, l[x] - 1, 1[ 

u(±l,y) = u(x,±l)=0 x,y€[-l,l] { ' 

with f(x,y) = 2, for x > 0, and f(x,y) = 1 for x < 0. The correct viscosity 
solution is 

{2(1 - |y|) x > 0, \y\ > 1 + lnx 
-2ln(x) x > 0, \y\ < 1 + lnx (66) 
x<0. 

We show in Figure [3] our results. In this case the convergence in the 
uniform norm fails. Convergence in the integral norm L 1 as proved in Section 
[4] is confirmed by Table [3} 

5.3 Test 3: Shape-from-Shading with discontinuous bright- 
ness 

The Shape-from-Shading problem consists in reconstructing the three di- 
mensional shape of a scene from the brightness variation (shading) in a 
greylevel photograph of that scene. The study of the Shape-from-Shading 
problem started in the 70s (see [T7] and references therein) and since then 
a huge number of papers have appeared on this subject. More recently, the 
mathematical community was interested in Shape-from-Shading since its for- 
mulation is based on a first order partial differential equation of Hamilton- 
Jacobi type (see [261 12S] ) . 
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Figure 3: Test 2. 



The equation related to this problem is the following: for a brightness 
function I(x,y) : M 2 D f2 — > [0, 1], in the case of vertical light source is 
vertical, to reconstruct the unknown surface, we need to solve 



\D U (x,y)\=(Jj^-lj, (x,y)eSl. (67) 

Points (x, y) where I is maximal (i.e. equal to 1) correspond to the particular 
situation when the light direction and n are parallel. These points are 



usually called "singular points" and, if they exist, equation (67) is said to 
be degenerate. The notion of singular points is strictly related to that of 
concave/convex ambiguity, we refer to |22| I20j for details on this point. 

It is important to note that, whatever the final equation is, in order to 
compute a solution we will have to impose some boundary conditions on dVt 
and/or inside O. A natural choice is to consider Dirichlet type boundary 
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Ax = h 


II ' 1 1 oo 






Ord(L x ) 


0.2 


1.0884 




0.4498 




0.1 


1.0469 




0.2444 


0.88 


0.05 


1.0242 




0.1270 


0.9444 


0.025 


1.0123 




0.0628 


0.9708 


0.0125 


1.0062 




0.0327 


0.9867 


0.00625 


1.0031 




0.0221 


0.5652 



Table 3: Test 2: experimental error. 



conditions in order to take into account at least two different possibilities. 
The first corresponds to the assumption that the surface is standing on a flat 
background, i.e. we set u(x, y) = for (x, y) £ dQ. The second possibility 
occurs when the height of the surface on the boundary (silhouette) is known: 
u(x, y) = g(x, y) for (x, y) E dCl. The above boundary conditions are widely 
used in the literature although they are often unrealistic since they assume 
a previous knowledge of the surface. 

Let us focus on two important points: 

• We note that a digital image is always a discontinuous datum. Is is a 
piecewise constant function with a fixed measure of his domain of reg- 
ularity (pixel) . So this is the interest of our analysis for discontinuous 
cases of /. 

• In the case of maximal gray tone (I(x) = 1), we are not in the Hypoth- 
esis introduced previously. In particular we have that / = in some 
points. We overcome this difficulty, as suggest in [ID]. We regularize 
the problem making a truncation of /. It is possible to show that this 
regularized problem goes to the maximal subsolution of the problem 
with e — > + . And that this particular solution is the correct one from 
the applicative point of view. 

We consider, now a test with a precise discontinuity on /, and we will 
discuss some issue about this case. 

We firstly consider a simple problem in ID to point out an aspect of the 
model. Let the function / be 

r Vl -x 2 if - 1 < x < 0.2 
1=1 ^ if 0.2 < x < 1 (68) 

I, otherwise 

we can see that we have a discontinuity on x = 0.2; despite this, because 
of the non degeneracy of the dynamics, the solution will be continuous. For 
this reason we can see that changing the boundary condition of the problem, 



Approximation of a discontinuous Eikonal Equation 



21 



the solution will be the maximal Lipschitz solution that verifies continuously 
the boundary condition. 




-1 - 



Figure 4: Sfs-data and solution with various boundary values. 



To see this we have solved this simple monodimensional problem with 
various Dirichlet condition, in particular we require u(—l) = 0, and u(l) = 
{-1,0.5,0,0.5,1}. With Ax = 0.01 and At = 0.002, we obtain the results 
shown in Figure |4} 

We can realize, in this way, an intrinsic limit of the model. It can not 
represent an object with discontinuities. We make another example that is 
more complicated and more close to a real application. 

We consider a simplified sfs-datum for the Basilica of Saint Paul Outside 
the Walls in Rome, as shown in Figure [5} We have not the correct boundary 
value on the silhouette of the image and on the discontinuities, so we impose 
simply u = on the boundary. Computing the equation with At = 0.001 
we get the solution described on Figure |6j 

We can see that, although the main features of the shape as the slope 
of the roofs, the points of maximum are well reconstructed. Despite it, the 
shape which we get is not so close to our expectations. We can try to get 
better results adding the correct height of the surface along the silhouette 
as discussed above and, in this case, we get the solution shown on Figure 
[7j We can notice a more convincing shape, but also in this case it is quite 
not satisfactory. For example we have that the correct boundary conditions 
we imposed are not attained, and we create some discontinuity on some 
parts of them. This is due to the fact that they can be not compatible 
with the statement of the problem. Essentially the limit which we can see, 
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Figure 5: Basilica of Saint Paul Outside the Walls: satellite image and 
simplified sfs-datum. 



as described above, is that we cannot have discontinuity on the viscosity 
solution (Theorem [2]) . 

We propose a different model for this problem, which allows discontinu- 
ous solutions. At this point we do not care about the physical interpretation 
of it, instead we are trying to find a solution closer to the correct solution. 
We want to solve the equation 

2 



max < —Du(x) ■ V a au(x,y) > = */-rn — \ — 1 x £ £1 ,„,,. 
\a\<i\ y ' fc =! V '} V I2 ( x <v) (69) 

u(x) = g(x) x S dQ 

with the map a : £1 — > M 2 ' 2 is 

(l + \I(x-h,y)-I(x + h,y)\)- p 

(l + \I(x,y-h)-I(x,y + h)\)- p 

(70) 

where p G K is a tuning parameter. Obviously this choice of the anisotropic 
evaluator a is a bit trivial. This pick is done for the sake of simplicity. More 
complicated proposal can be found for example in pQ. 

In this way we use the results about the degeneracy of the dynamics 
permitting to the viscosity solution to be discontinuous. Of course this is, 
in some sense, the opposite situation with respect to the classical formula- 
tion: in this case every non smooth point of the surface is interpreted as 
discontinuity and we try to reconstruct it using the data coming from the 
silhouette. 

The results are shown in Figure [8] and in Table [4] we can see an accuracy 
comparison of the various procedure. 
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. 4 




Figure 6: Test 3: reconstructed shape without boundary data. 



Test 


II ' 1 loo 


• L 1 


w/o correct boundary data 


1.7831 


1.5818 


w boundary data 


0.8705 


0.5617 


w boundary + disc, detect. 


0.7901 


0.3062 



Table 4: Test 3: Comparison between various methods 
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Figure 7: Test 3: Dirichlet condition on the silhouette. 




Figure 8: Test 3: Dirichlet condition and discontinuous dynamics. 



